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Abstract. We present a straightforward and reliable continuous method for 
computing the full or a partial Lyapunov spectrum associated with a dynamical system 
specified by a set of differential equations. We do this by introducing a stability 
parameter (3 > and augmenting the dynamical system with an orthonormal k- 
dimensional frame and a Lyapunov vector such that the frame is continuously Gram- 
Schmidt orthonormalized and at most linear growth of the dynamical variables is 
involved. We prove that the method is strongly stable when (3 > — where A^ is 
the fc'th Lyapunov exponent in descending order and we show through examples how 
the method is implemented. It extends many previous results. 
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1. Introduction 

One of the most striking coordinate independent characterizations of a compact 
dynamical system of dimension d is its Lyapunov spectrum. It associates to each 
orbit of the system a set of d real values which describes exponential instabilities of 
infinitesimal deviations from the orbit. Furthermore, for an ergodic dynamical system 
this set almost surely (in the measure theoretical sense) does not depend on which 
orbit you consider. More precisely, given a smooth vector field v in R d we look at 
the d-dimensional evolution equation x = v(x). For an initial condition x(0) = xq we 
integrate this to obtain a corresponding orbit x(t) = 0*(xo). The stability of such an 
orbit can then be examined by looking at the evolution of a nearby orbit x(t) +u(t) and 
linearizing the equations of motion in u : 

*(*) = = j(x(t)) u (t). (i) 

Integrating this along the orbit we obtain the tangent map u(t) = M Xo (t)uo in which 
the transition matrix M Xo (t) = d(p t (xo)/dxo is a d by d matrix valued function of time. 
The exponential instabilities of a trajectory are now reflected in its eigenvalue spectrum 
or rather the spectrum of the symmetric product 

Ml(t)M XQ (t). (2) 

The spectrum of this matrix is real and positive and we order it as follows : 

til(t) > t4(t) > ... > > o. (3) 

Although the values a priori depends on the initial point xo chosen, one has : 

Theorem (Oseledec, [III 0]) : If ^ is an ergodic probability measure for the dynamical 
system then for ^-almost every xq : 

A fc = lim - log /i k (t) (4) 
exists and is independent of the initial point. 

In other words taking an arbitrary (with respect to an ergodic measure) initial point 
and calculating the above limits, with probability one you will get its unique Lyapunov 
spectrum, {Ai > A2 > ... > A^}. From its definition it is not difficult to show that 
the spectrum is independent of the choice of coordinate system and thus being intrinsic 
invariants, their existence and the above theorem is of fundamental importance in the 
theory of dynamical systems. Determining the spectrum (or part of it) has in fact grown 
into a small industry in modern non-linear physics. 

From the practical or numerical point of view the above description is insufficient 
as the matrix M T (t)M(t) pretty fast becomes singular since its eigenvalues separate 



3 



exponentially in time (assuming that not all Lyapunov exponents are equal) thus making 
it difficult to measure the spectrum. Now, several methods have been developed in 
order to overcome this problem. Fr0yland [[J uses a systematic but rather complicated 
evolution equation for co-matrices which still has exponential growth, although in a 
controllable way. Meyer M makes use of symplectic transformations to derive the 
spectrum for Hamiltonian systems also involving exponential growth. Habib et al. [0] 
uses a particular representation of the Lie algebra so(2) to obtain an evolution equation 
for a Hamiltonian system in 2 dimensions (which appears to work also in dimension 4) 
which only involves linear growth but which is rather complicated and depends on the 
Hamiltonian form of the dynamical system. Bennetin et al. || (for one exponent) and 
later Shimada et al. |7j (for the whole spectrum) suggest to renormalize (Gram-Schmidt) 
at regular intervals of time a set of stability vectors picking up the exponents during 
the renormalization procedure. This method works well and is often used in practice 
when calculating the spectrum (an example is Gong || in which the authors considers 
the Lyapunov spectrum for a compact lattice Yang- Mills SU(2) theory). Goldhirsch et 
al. || present a continuous version of this procedure (cf. below) and they develop a 
set of differential equations for the eigenvalues and eigenvectors of the stability matrix 
M Xo (t) itself, a method, however, unsuited in presence of a degeneracy of eigenvalues. 

Here we shall present a unified approach in which we augment the dynamical system 
with an orthonormal frame and a Lyapunov vector such that the augmented system is 
dynamically strongly stable and involves at most linear growth and such that the Lya- 
punov spectrum is obtained almost surely (in the measure sense of choosing an arbitrary 
initial point and frame). The method is not constrained to Hamiltonian systems, it ap- 
plies to any finite dimensional dynamical system and is straightforward to implement 
on a computer. 

2. Continuous Gram-Schmidt orthonormalization 

We define a time dependent orthonormal fc-frame to be a set of k (k < d) orthonormal 
vectors : 

£ (t) = {ei(t), ..,e k (t)}, (ei,ej) = Sij, l<i,j<k, (5) 

where (•, •) is the usual Euclidian product in R d . Using this frame we let Ji m = (ei, Je m ) 
denote the matrix elements of the Jacobian matrix J and we note that these ma- 
trix elements depend on time both through the Jacobian and the frame. We intro- 
duce a stability parameter (3 > and define the (symmetric) stabilized matrix ele- 
ments L mm = J mm + (3((e m ,e m ) - 1) and L im = J lm + J ml + 2/3(e h e m ). Finally, let 
A = {Ai(t), Afc(t)} be a fc-dimensional real vector. 
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The augmented dynamical system is now given by the following set of differential 
equations (of which the first two are vector equations) : 

x = v(x) , 

e m = Je m -J2i<m e iLi m m = l,...,k, (6) 

A m J mm 1, /c . 

For the dynamical evolution of these equations we have : 

Theorem : Let xq be an initial point for which the associated Lyapunov spectrum ( cf. 
equation^) X% > \ 2 > ... > A^ exists. Set A(t = 0) = 0. Choosing the stability parameter 
(3 > — Afc then for almost any (i.e. with probability 1 when choosing randomly) initial 
frame £(t = 0) the time evolution of the dynamical system (Qj yields : 

Km ~A m (*) = Am, m = 1, k. (7) 

t— »oo t 

Thus, by following a trajectory of the augmented system we obtain almost surely the 
k first elements in the Lyapunov spectrum for the given orbit. The somewhat peculiar 
condition on the stability parameter is satisfied e.g. by setting /3 > max|| e || =1 ( — (e, Je)) 
where the maximum is over all unit-length vectors e and over the relevant region of phase 
space. Dynamically such a choice corresponds to finding the strongest local contraction. 

The proof of the Theorem is given in appendix A where in particular, it is shown 
that the dynamics preserve orthonormality of the frame. When the elements of J are 
assumed bounded in phase space we see that the above dynamical system only involves 
at most linear growth of the dynamical variables (through the A m 's). Given a dynamical 
system with an ergodic measure we see by combining the Oseledec Theorem and the 
above that the Lyapunov spectrum is obtained almost surely by choosing an arbitrary 
initial point and an arbitrary initial frame. 

An interesting case is when k = d, i.e. when we want to calculate the complete 
spectrum. In this case, our orthonormal frame is complete and we may expand Je m in 
equation @ in terms of the basis vectors themselves. Setting (3 = we get 

C m = J! e lJlm — X! e lJrnl = e l^-lm: (8) 
l>m l<m I 

where Ai m is anti-symmetric and thus by construction a generator of orthonormal trans- 
formations (of our frame). A straight-forward linear analysis shows that the resulting 
dynamical system is marginally stable. 

3. Numerical results 



In the following we apply the above method in order to calculate the Lyapunov spectrum 
for two standard systems; the Lorentz system and a 3 degrees of freedom Hamiltonian 
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system with a quartic potential. In order to get good statistics on the Lyapunov 
exponents we determine a time over which to integrate the systems in order to get 
reasonable convergence and then make 1000 runs of each system with random initial 
conditions. 



The Lorentz system [ FD | is given by: 
x = — ax + ay 

y = rx — y — xz (9) 
z = xy — bz 

with the usual choice of parameter a = 10, r = 28 and b = 8/3. In figure |]we present 
the resulting Lyapunov exponents of a single run for the Lorentz system. The result of 
1000 runs is given in table |1]. 



Table 1. The average of the Lyapunov exponents of 1000 runs of the Lorentz system 
and their root mean square deviations. The sum of the exponents is —13.6667 = 
—a — 1 — b as expected from equation (9). 
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rms dev. 
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1.4-10" 5 


8.3-10- 4 
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-14.5724 


4.6-10- 3 



Next we apply our method to a Hamiltonian system with a quartic potential: 

H = Pl+P 2 y+P 2 z x 2 y 2 + y 2 z 2 + z 2 x 2 x 4 + y 4 + z 4 

2 2 32 [ ' 

The last term is added in order to have a compact phase space and avoid having to deal 
with problems of convergence related to near-integrable motion along the coordinate 
axes, but still chosen with a sufficiently small prefactor in order not to stabilize the 
dynamics. A single run of the resulting 6-d system is presented in figure |2| and the 
corresponding averages of 1000 runs are given in table 

Looking at the results for the two systems there is one striking difference which 
is apparent both for the single run results shown in the figures and in the averages 
given in the tables. Whereas the (finite-time) exponent A2 for the Lorentz system, 
corresponding to the marginally stable direction along the flow, fluctuates around zero, 
the equivalent (finite-time) exponent, A3, for the Hamiltonian system is clearly different 
from zero though converging for increasing t. The difference can be made even more 
apparent when plotting instead e A * as in figure |3|. e A * is the stability eigenvalue for 
the marginally stable direction over the entire integration. For the Lorentz system 
this remains constantly close to one as expected by a marginal eigenvalue. For the 




-14.5 



-14.54 



-14.58 



100 200 300 400 500 600 700 800 900 1000 

t 




100 200 300 400 500 600 700 800 900 1000 

t 



Figure 1. The 3 (finite-time) Lyapunov exponents from a single run of the Lorentz 
system. 



Table 2. The average of the Lyapunov exponents of 1000 runs of the Hamiltonian 
system (^1) and their root mean square deviations. 
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Figure 2. The 3 positive (finite-time) Lyapunov exponents from a single run of (JToj) . 
We present only the positive exponents since, due to the symplectic structure of the 
equations, the absolute value of the negative exponents are exactly equal to the positive 
ones. 



Hamiltonian system, however, it grows linearly. The explanation lies in the degeneracy of 
this eigenvalue. The Hamiltonian (|T0|) has another marginally stable direction associated 
with gradif . In local coordinates the Jacobian for the flow will in general take on a 
Jordan normal form for the two marginally stable directions, i.e. we can not expect 
to find a full set of eigenvectors for the Jacobian. One may consider the possibility of 
factoring out the marginally stable directions and not include them in the integration, 
but here we just note that we know that the corresponding exponents exactly equals 
zero for infinite time. 
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Figure 3. e where A is the exponent corresponding to the marginally stable direction 
along the flow for (a) the Lorentz system and (b) the Hamiltonian system ([h]) . 

To illustrate the effect of the stabilizing factor (3 we consider the time to reach a 
certain error level on the orthonormality of the basis S with varying (3. To be more 
specific we consider the error e = [J2ij=i,k(( e i^ e j) ~ <%) 2 ] an d the time T(0) it takes 
for this error to grow to a certain level as a function of (3. Since this will depend on 
the initial conditions we have taken the average of 10 runs for each of the systems. 
The result is presented in figure [|. We let the systems run a maximum time before 
deciding that the chosen error level would not be reached. For the Lorentz system this 
time was set at 2000 and for the Hamiltonian system at 10000. The curves of figure 
[| therefore saturate at these values. For the Lorentz system T({3) increases drastically 
near (3 = —A3 as expected, whereas for the Hamiltonian system the picture is a little less 
clear with the large increase of T{(3) happening at a slightly larger value than f3 = —Xq. 
The difference is due to the relative high homogeneity of the Lorentz system vs. the 
rather strong dependence of the local finite-time stability exponents on the phase space 
position in the Hamiltonian system. Based on these results we have chosen to set /? = 20 



9 



for the Lorentz system and (5 — .5 for the Hamiltonian system. The point here is to set 
(3 sufficiently high to stabilize the given system under integration, but not excessively 
high since this could easily lead to unnecessarily high requirements on the accuracy of 
the integration routine. 

The method to calculate Lyapunov exponents of ODE's that we have presented 
in this paper is nothing but a continuous version of the standard Gram-Schmidt 
orthonormalization procedure. As mentioned above this was already proposed by 
Goldhirsch et al. [§, eqs. (5.12) and (5.26)] but without the crucial stability term. 
Apart from the aesthetic pleasure of formulating the whole procedure as one set of 
differential equations, the method will show its usefulness when calculating Lyapunov 
spectra where the difference between the largest and the smallest exponent is large. 
In such a case, using standard methods, one rapidly loses accuracy on the eigenvector 
associated with the lowest exponent and therefore also of the exponent itself. One 
would therefore have to employ the Gram-Schmidt orthonormalization quite often; the 
continuous orthonormalization naturally avoids this problem with the /3-term replacing 
the stabilizing effect of the (non-linear) Gram-Schmidt procedure. On the other hand, 
if one is only interested in calculating the largest exponent for a given system there 
is essentially no difference between standard methods and the method given here 
since orthonormalization is unnecessary except to avoid a possible numerical overflow. 
Computationally the continuous method is somewhat heavier than standard methods. 
To compute k exponents of a d- dimensional system one needs, in addition to the usual 
Je m operation (0(kd 2 )), to compute eiJe m and e«e m , both 0(k 2 d). For a full spectrum 
the computation is thus a factor of 3 heavier, whereas for partial spectra it will be 
somewhat less. 
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Figure 4. The time T it takes to reach an error level of 10 3 on the orthonormality 
of the basis £ for the Lorentz system (a) and the Hamiltonian system (p"o|). 



Appendix A. 

Appendix A. 1. Proof of the Theorem 

We first remark that when the frame is orthonormal the stabilizing terms vanish 
identically. We shall show that the resulting equations correspond to a differential 
version of a Gram-Schmidt orthonormalization of a set of k independent tangent vectors 
evolving in time according to equation (p]). 

Thus, consider a time dependent set of vectors {fi,---,fk} which are linearly 
independent and satisfies f m = J(x(t))f m . We expand this set of vectors uniquely 
into an orthonormal set {ei, e&}, i.e. f m = J2i<m e i K im where {ni m }i< m are a set of 
time dependent parameters with positive diagonal elements, i.e. K mm > 0. 
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We shall show that the e m 's and the diagonal elements K mm satisfy the differential 
equations : 

Cm J&m &mJmm Si<m &l\Jlm J'ml) i / All 



' mm " in in 



We prove it by induction. So assume it is true for the vectors {e 1; ...,e m _i}. From 
fm = Jfm we get : 

Y {zi K lm + e l^lm - JeiKi m ) = . (A2) 
l<m 

By orthonormality (e m , e m ) = and by the induction hypothesis (e m , e/) = J m i for 
I < m so we obtain by taking the scalar product with e m : 

^""^ {Jml^lm) 4" K mm Jml^lm = (A3) 

Z<m l<m 

and hence that K mm = J mm n mm . Again for I < m by the induction hypothesis ei — Jei 
and hence also J2i<m(^i K im+€-iki m — JeiKi m ) are in the span of ex, e m -i (just re-arrange 
the terms in (|A1[)) and hence, we may rewrite equation (IAI3) as follows : 



^m^mm ^-mJmm^mm J&m^mm ^ t j (A4) 

Z<m 

where q are some (time dependent) parameters. For I < m the scalar product (e;,e m ) 
is constant (zero) in time. Whence (e;,e m ) = — (e m ,e;) = — J m ; and we may take the 
scalar product with t\ in (|A4|) to obtain 



{Jml "I" Jlm)K mm — C\ . (A5) 



Finally, inserting this in (|AJ) and dividing by K mm we get the desired equation for e r 
From (pi) and the above we get the relationship K mm (t) = exp(A m (t))ft mm (0). 



Appendix A. 2. The Lyapunov spectrum 

Next, we will show that provided the Lyapunov spectrum exists (we assume so from 
now on), the limits lim^oo jA m (t) will almost surely give the spectrum (in descending 
order). In order to do this we take the positive matrix M T M from equation (0) and 
diagonalize it to obtain : 

M T M = Y, &rfl™ ® a m (A6) 

m 

where {ak} m =i,..,k is a set of orthonormal vectors and the /z^'s are as in ([]). One has 
||Mw|| 2 = f I k( a ki u ) 2 which geometrically means that the image of a sphere ||w|| = 1 
will be an ellipsoidal with principal axes {/!/«}. Now let 1 > r > be a fixed constant 
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and let u be a unit length vector. Suppose that \(ai, u)\ > r > at all times. Then we 
have : 

l4 > || Mu || 2 > r 2 /i 2 (A7) 
and hence it follows that : 

lim ^log||Mu|| 2 -Ai= Urn - logAii-Ai + 7(log||Mu||-logAii) = 0(A8) 

since Ai = lim | log /xi and the expression in the parenthesis is uniformly bounded : 
log(r) < log ||Mm|| — log/ii < 0. Now, the vector a\ actually depends on time, but since 
u is chosen at random we have at any given instant that |(oi, u)\ > r > with a prob- 
ability p(r) which tends to 1 as r tends to zero. This follows from simple geometrical 
considerations on the area of the d-ball, compared to the part of it for which the above 
inequality holds. Hence the above results holds with probability p(r) and since r > 
was arbitrary it follows that with probability 1 the limit of 1/t log ||Mm|| will exist and 
be the maximal Lyapunov exponent. Inserting here u = ei(0) and Mu = exp(Ai(t))ei(t) 
we obtain the desired result. 



The general case is shown by considering growth rates of exterior products. We 
shall only show the explicit calculations for the first two Lyapunov exponents, noting 
that the formulae easily generalize. 

We recall that if {e 1: ...,e d } is a basis for V, then the formal exterior products 
{ek A ei} k< i is a basis for the vector space A 2 V = V A V. One may define the scalar 
product : 



det 



(ui,vi) {u u v 2 ) 

{U2,vi) (U 2 ,V 2 ) 



(ui A u 2 ,vi A v 2 ) 

as well as the action of A 2 A = A A A : 

A A A(ui A u 2 ) = (Aui) A (Au 2 ) . 

Consider now the action of the matrix A 2 (M T M) in the following way : 

(it! A u 2 , A 2 (M T M)v 1 A v 2 ) = 

(u! A u 2 , (M T Mv 1 ) A (M T Mv 2 )) = 

det{(w;, M T Mvj)} = 

det {(M u h Mvj)} = 

(A 2 Mu x A u 2 , A 2 Mv x A v 2 ) = 

(ui A u 2 , (A 2 M) T (A 2 M)v 1 A v 2 ). 



(A9) 



(A10) 



(All) 



In particular, this shows that A 2 (M T M) = (A 2 M) T (A 2 M), an identity which will allow 
us to estimate the growth rate of the product of the two largest eigenvalues of M. Using 
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the diagonalization (|A6| ) as well as linearity and anti-symmetry of the wedge product 
we get : 

A 2 (M T M)(u 1 A u 2 ) = ^f^f^cii A ajfa A a^U\ A u 2 ) . (A12) 

i<j 

In particular : 

(«iA« 2 ,A 2 (M t M)(«iAm 2 )) = Y,tffj%(a i Aaj,u 1 Au2) 2 ■ (A13) 

i<j 

and using (|Allj) we deduce the inequality : 

A*iA*a| (oi A a 2 , «i Am 2 )| < || A 2 M(u\ A u 2 ) || < /ii£t 2 • (A14) 
We can then repeat the arguments from above to show that with probability 1 : 

lim - log || A 2 M(ui A u 2 ) || = Ai + A 2 (A15) 

and by anti-symmetry of the wedge the left hand side will apart from a uniformly 
bounded contribution (which almost surely vanishes in the limit) equal 

lim-log||/iA/ 2 || =lim-log|«ai(t)K2a(*)| = lim-(A 1 (t)+A 2 (t)), (A16) 

and using the already obtained formula for the exponent Ai the desired result follows 
for A 2 . 



Appendix A. 3. Linear stability theory 

The results obtained above are numerically reliable provided the frame stays 
orthonormal during the time evolution. In particular the variables A [m = (ej, e m ) — 5i m , 
1 < I, m < k should all vanish. By straight-forward differentiation one verifies that they 
satisfy the following set of differential equations : 

A pm = —{2(3 + L mm + L pp ) — ^2 ApiLi m — Ai m Li p . (A17) 

l<m l<p 

It is clear that A pm = 0, for all p and m, is a fixed point of these equations. In order to 
analyze its stability we substitute A pm — > A pm + 5 pm and linearize in the variables S pm 
to find : 

Spm = ~ (2/3 + Jmm + Jpp)b~pm + C({(5pj}/< m , {5/ m }/< p ) (A18) 

where G is linear in the 5's but depends only on the preceding variables. Here we have 
used the natural lexicographic ordering : (11) < (21) = (12) < (22) < (31) = (13) < 
(32) = (23) < ... It follows that the stability of these equations is determined only by 
the stability of the first term, i.e. of the differential equations (for all p and m) : 

z = -(2(3 + J mm + Jpp)z. (A19) 
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This equation happens to be analytically solvable, surprisingly in terms of our Lyapunov 
vectors themselves, 

z(t) = exp(-2/3t - A rnm (t) - A PP (t)), (A20) 

and we see that stability is ensured provided 

P > - lim -A mm (t) = -X m , (A21) 

t~+co t 

for all m = l,...,k. As our A's are ordered decreasingly it suffices to have (3 > — 
From equation ( |A19| ) we also see that stability follows by using the more conservative 
bound obtained by setting (3 > max|| e || =1 (— (e, Je)). 
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